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I.  Introduction 


.  1 


Computer  simulation  of  condensed  phases  and  complex  phenomena 
adds  a  new  important  dimension  to  physics  research.  In  particular,  it  has  proven 
invaluable  in  bridging  the  gap  between  experiment  and  theory  in  those  cases 
where  either  (or  both)  of  the  two  is  (are)  inaccurate,  insufficient  or 
fundamentally  limited.  Moreover,  computer  simulation  vastly  extends  the  scope 
of  model  building  and  testing  from  the  traditional  domain  of  linear  equations 
and  analysis  to  that  of  nonlinearity,  which  is  more  a  rule  than  an  exception  in 
nature. 

There  are  actually  a  multitude  of  numerical  techniques  often  vaguely 
referred  to  as  computer  simulation  or  computer  experiments.  They  range  from 
numerical  solutions  of  partial  differential  equations  to  Monte  Carlo,  Langevin 
and  molecular  dynamics  simulations  to  path  integral  evaluations.  Likewise,  the 
physical  applications  include  such  diverse  systems  as  turbulence  and 
hydrodynamics,  lattice  gauge  theories  of  quantum  chromodynamics,  phase 
transitions  and  critical  phenomena,  stochastic  dynamics  of  growth  and  pattern 
formation,  and  electronic,  structural  and  thermodynamic  properties  of  materials. 
This  talk  will  concentrate  on  the  last  one.  There  has  been  imoressive  progress’ 
in  that  area  during  the  last  few  years,  both  in  the  development  of  algorithms  and 
computer  codes  and  their  implementation  to  large-scale  computing,  as  well  as  in 
the  development  and  refinement  of  useful  phvsical  and  mathematical  concepts 
and  devices,  in  particular,  molecular  dynamics  simulation*-’  combined  with  a 


proper  description  of  electron-mediated  forces,  provides  a  powerful  technique 
for  a  computer-based  microscopic  description  of  interesting  and  important 
material  properties.  I  shall  discuss  (i)  the  recent  ideas  on  how  to 
systematically  go  beyond  the  rigid  (pair)  potential  description  of  interatomic 
interactions  and  make  contact  with  electronic  structure  calculations  for 
low-symmetry  systems  and  (ii)  present  selected  applications  of  molecular 
dynamics  simulations  to  physical  problems  of  current  interest. 

2.  Beyond  classical  potentials 

Molecular  dynamics  simulations  are  conceptually  simple  yet 
increasingly  useful  in  elucidating  the  properties  of  matter  in  the  microscopic  to 
mesoscopic  scale.  The  usefulness  and  quantitative  reliability  ,  nowever,  often 
hinge/  on  to  what  extent  the  force  laws  between  interacting  atoms  and 
molecules  can  be  specified.  While  relatively  simple  force  laws  are  feasible  for 
such  systems  as  rare  gases  and  polar  molecules,  descriptions  based  on 
few-body  forces  are  not  universally  applicable  in  metallic  and  covalently  bonded 
matter.  (Such  descriptions  may,  however,  be  useful  in  limited  applications  3S 
evidenced  by  the  vast  body  of  work  on  pair  potential  simulations  for  metals^  and 
the  recent  efforts4  in  modelling  Si  and  Ge  Dy  two-body  and  three-body  forces  ) 

Here  I  discuss  two  recent  approaches  to  interatomic  interactions  The  r  >rst  is  a 
scheme  to  perform  first-principles  electronic  structure  calculations  in  parallel 


•3. 


to  the  atomic  dynamics,  while  the  second  entails  a  class  of  more  approximate 
but  computationally  fast  methods  of  obtaining  the  force  laws. 

2.1.  Simulated  annealing  and  related  methods 

A  few  years  ago  algorithms  based  on  thermodynamic  concepts  were 
introduced  to  solve  large-dimensional,  nonlinear  optimisation  problems,  in  a 
short  period  of  time  these  techniques,  the  most  famous  among  which  is  the 
simulated  annealing5  method,  have  attracted  a  good  deal  of  attention 
Optimisation  problems  involve  a  "cost  function",  which  has  to  be  minimized  with 
respect  to  the  relevant  parameters  (degrees  of  freedom).  Simulated  annealing  is 
cased  on  the  analogy  between  minimizing  the  cost  function  in  an  optimisation 
problem  and  the  physical  problem  of  obtaining  a  ground  state,  e  g  the  process  of 
growing  a  good  crystal  by  slow  cooling  from  its  melt.  In  order  to  avoid  defects  or 
glassy  areas  in  the  crystal,  slow  annealing  is  necessary.  This  prohibits  the 
system  from  getting  stuck  into  a  local  free  energy  minimum  In  the  general 
optimisation  problem,  this  may  be  accomplished  by  using  eg.  the  Metropolis 
Monte  Carlo  algorithm  to  accept  random  displacements  in  the  parameter  space 
The  temperature  enters  as  a  control  parameter  in  the  Boltzmann  factor,  wmcn  is 
adjusted  to  ensure  the  passage  to  the  global  minimum  As  a  last  step,  a 
steepest-descent  quench  can  be  used  to  lock  into  the  global  minimum 

A  very  important  contribution  to  the  simulated  annealing  ideas  for 
electronic  ana  structural  properties  is  due  to  Car  and  Parnneilo^  They  compined 


1C 


the  electronic  density-functional  Hamiltonian  with  the  ionic  kinetic  energy  into 
a  cost  function  to  be  globally  minimized.  Traditionally,  in  electronic  structure 
calculations  one  first  minimizes  the  electronic  energy  for  fixed  ionic  positions. 
According  to  the  Born-Oppenheimer  principle,  the  (Hellmann-Feynman)  forces  on 
ions  can  then  be  evaluated  and  the  system  is  allowed  to  relax.  This  leads  to  a 
series  of  elaborate  self-consistent  calculations  with  the  necessary  matrix 
diagonalizations.  Even  with  the  biggest  available  computers,  such  a  procedure 
becomes  orombitively  expensive  for  systems  containing  more  than  a  dozen  or  so 
atoms,  as  the  number  of  degrees  of  freedom  to  be  optimised  increases  this  is 
so  despite  the  great  simplifications  provided  by  3t>  initio  pseudopotential  theory 
and  local-density  approximations  for  exchange  and  correlations  Another 
bottleneck  in  these  calculations  is  the  cost  of  large  matrix  diagonalisations, 
which  scale  as  proportional  to  M~\  where  M  is  the  total  number  of  basis 
functions  needed  to  represent  the  electronic  states  and  charge  densities 

Car  and  Parnnello  suggested  that  the  combined  cost  function  of 
electrons  and  ions  be  minimized  by  using  simulated  anneling  in  the  phase  space 
of  nuclear  coordinates  and  density-functional  single-particle  electronic 
wavefunctions  Rather  than  doing  a  Monte  Carlo-type  search,  they  add  a  fictitious 
classical  noetic  energy  '.with  an  arbitrary  mass)  to  the  electronic  fields,  and 
derive  the  lagrangian  equations  of  motion  for  both  electronic  and  ionic  degrees 
of  freedom  Additional  dynamics  can  be  assigned  to  eg.  the  volume  and  snape  of 
tne  system  Orthonormality  of  the  electron  states  is  imposed  by  the  technique  of 
Lagrange  multipliers  Tne  result  is  a  set  of  coupled  equations  of  motion,  which 


can  be  simultaneously  solved  for  the  ion  coordinates  and  the  electron  fields 
using  standard  discrete  time-step  methods  of  classical  molecular  dynamics. 
Starting  from  a  suitable  set  of  initial  conditions,  the  "temperature"  associated 
with  the  fictitious  classical  dynamics  of  electrons  can  be  gradually  lowered  so 
that  the  system  can  relax  to  its  global  ground  state.  By  adjusting  the  electron 
classical  "mass'  one  can  (and  should)  stay  in  the  vicinity  of  the 
Born-Oppenneimer  surface. 

This  technique  has  two  important  merits.  Firstly,  it  provides  a 
systematic  and  practical  way  of  combining  real  molecular  dynamics  of  the  ions 
with  interactions  truly  derived  from  electronic  structure,  i.e.  no  assumptions  of 
the  force  laws  are  required.  This  is  very  important  for  looking  at  such  delicate 
things  as  structural  phase  transitions,  gram  boundary  structures7,  finite 
temperature  structures  of  clusters8  etc.  Secondly,  as  the  simulated  annealing 
technique  expicitly  avoids  direct  matrix  diagonalisation  steps,  its 
computational  cost  for  the  electronic  part  of  the  problem  increase  only 
proportionally  to  MN  (N  is  the  number  of  electrons)  as  far  as  the  "iterative 
diagonalisation"  is  concerned.  There  are  other  parts  of  the  calculation  such  as 
the  orthonormalisation  of  the  occupied  states  (MN^),  the  generation  of  the  cnarge 
densities  (MN),  and  the  solution  of  the  Poisson  equation,  which  scales  as  MlnM 
(using  fast  Fourier  transforms)  However,  it  is  obvious  that  the  Car-Parrinello 
technique  will  be  the  one  of  choice  for  very  large  electronic  structure 
calculations. 

in  fact,  if  one  freezes  the  nuclear  positions  and  only  looks  at  the 


electronic  minimization  problem,  there  probably  are  even  more  efficient  dynamic 
method  for  the  purpose  of  iterative  diagonalisation.  The  cost  function  is  in  this 
case  convex  (at  least  for  nonmagnetic  systems;  magnetic  ones  can  exnipit 
bistability)  and  steepest-descent  techniques0  are  feasible  Another  recent 
suggestion10  deais  with  speeding  up  the  orthonormalisation  step  By  using 
Schrodinger  dynamics  and  a  Fourier  bandpass  filter  one  can  quickly  sieve  out  the 
eigenfrequencies  and  associated  eigenfunctions.  For  evaluating  the  total  electron 
density,  which  usually  is  the  quantity  of  interest  in  self-consistent 
calculations,  one  can  show  that  the  procedure  can  be  optimised  to  scale  as 
MN3/2  . 

Another  aspect  of  the  computational  cost  of  iarge-scale  electronic 
structure  calculations  concerns  the  basis  sets  used  to  represent  the 
wavefunctions  and  densities.  Two  commonly  used  basis  sets  in  large  calculations 
are  plane  waves  and  Gaussians.  The  latter  do  not  constitute  a  complete  and 
orthonormal  set,  and  consequently  quadruple  sums  over  products  of  the  basis 
functions  are  needed  in  evaluating  some  of  the  energy  terms.  This  quickly 
becomes  expensive.  Using  plane  waves  is  equivalent  to  using  an  equally  spaced 
grid  in  three  space  dimensions,  if  one  needs  a  resolution  good  enough  to  describe 
core  electrons  in  atoms,  millions  of  plane  waves  are  needed  to  represent  the 
simplest  unit  cells.  To  overcome  this  problem,  pseudopotentiais  can  be 
constructed  to  eliminate  the  core  orbitals,  if  these  potentials  are  local,  the 
computational  task  remains  relatively  manageable.  However,  to  describe  atoms 
sucn  as  oxygen  nonlocal  pseudopotentiais  are  necessary  This  leads  to  mixing  of 


the  Fourier  components  of  the  basis  with  a  considerable  increase  in  tne 
computational  difficulty. 

The  finite  element  method11  provides  in  principle  an  attractive 
alternative  as  a  basis  set.  In  this  method,  low-order  polynomial  shape  function* 
are  defined  on  a  real  space  grid  so  that  the  functions  are  strictly  zero  outside 
their  respective  elements.  The  mesh  can  be  varied  so  as  to  give  high  resolution 
where  needed.  If  core  orbitals  are  kept,  this  amounts  to  stretching  the  resolution 
near  nuclei  so  that  the  important  cusp  condition  can  be  fulfilled.  It  is  also  very 
easy  to  represent  any  function  in  terms  of  these  shape  functions.  The  method 
leads  to  a  very  sparse  (but  large)  Hamiltonian  and  overlap  matrices,  whicn  are 
nearly  optimally  amenable  to  iterative  schemes  such  as  the  Car-Parnnelio 
method.  Moreover,  fast  algorithms  (e.g.  the  multigrid  technique12)  are  avaiiapie 
for  solving  the  Poisson  equation  using  finite  elements.  Combining  the  methods 
listed  above  may  eventually  lead  to  breaking  the  bottleneck  in  very  large 
electronic  structure  calculations  with  real  atomic  and  molecular  dynamics 

2.2.  Effective-medium  theories 

Even  with  the  increase  in  computing  power  and  the  ingenious  new 
algorithms,  the  task  of  reliably  calculating  total  energies  for  interacting  atoms 
will  long  remain  problematic  in  eg.  studying  the  dynamics  of  complex, 
low-symmetry  extended  defects.  Another  limitation  of  many  nrst-pr-nciples 
calculations  is  that  the  amount  of  physical  insight  they  directly  provic  .»  often 


limited  In  order  to  identify  the  important  parameters  and  pnysical  trends 
extendable  to  other  systems  simpler  models  have  to  be  extracted  and  deveiooed. 

An  alternative  approach  to  calculating  the  total  energy  of  a 
complicated  system  is  provided  by  the  effective-medium-theory1  ^  in  its  various 
disguises.  The  simplest  definition  for  these  approaches  is  that  they  aim  to  go 
oeyond  the  pair-potential  description  of  condensed  matter  to  include 
density-dependent  interactions.  The  pair-potential  model  has  serious  drawbacks 
m  applications  to  energetics  of  metallic  and  covalent  solids  (  it  may  still  be 
very  useful  for  generating  qualitative  information  such  as  the  local 
configurations  of  atoms).  For  example,  unless  the  elastic  constants  satisfy  the 

Cauchy  relation  Cj2  =  ^44,  which  is  seldom  true  in  real  solids,  an  equilibrium 

pair  potential  cannot  reproduce  them.  One  suggestion  to  deal  with  this  problem  is 
to  add  to  the  total  energy  3  term  which  depends  on  the  macroscopic  volume,  by 
analogy  with  simple-metal  pseudopotential  perturbation  theory.  The 
volume-dependent  term  provides  a  fictitious  external  pressure,  which  balances 
‘.he  Cauchy  pressure.  It  is  however  clear  that  the  bulk  modulus  is  no  more 
uniquely  defined  (unless  the  volume-dependent  term  is  linear  in  volume).  This 
difficulty  may  lead  to  a  gross  mispresentatlon  of  cohesive  properties  so  that 
even  qualitative  simulating  of  such  processes  as  crack  opening  or  vacancy 
clustering,  which  involve  volume  changes,  become  suspect. 

The  basic  idea  behind  the  effective-medium  theories  is  simple:  the 
total  energy  of  any  given  atom  is  determined  by  the  effect  of  the  surrounding 
atoms  'n  the  spirit  of  density-functional  theory,  the  electron  density  r.(r)  is 


singled  out  as  the  key  variable.  One  writes  a  generically 
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Etot(N)-.2Ft(nt(Ri»  *  E6s  .  (1) 

The  first  term  is  the  "embedding"  energy  of  an  atom  to  its  (instantaneous) 
environment;  its  dependence  is  through  the  electron  density  of  all  the  other 
atoms.  An  ansatz  for  this  is 

nf<R)  =  2na(Rj  -R«),  (2) 

h*  ‘ 

where  na(R)  is  a  density  cloud  associated  with  an  atom  at  R.  The  first  term 

contains  the  leading  many-atorn  contributions  arising  from  the  interaction  of  the 
atom  with  the  electron  density  tails  of  all  the  other  atoms.  The  function  F  can  De 
calculated  for  a  suitably  chosen  reference  system.  A  particularly  convenient  one 
is  the  homogeneous  electron  gas,  for  which  the  embedding  energies.,  screening 
clouds  and  associated  electrostatic  potentials  can  be  calculated  exactly  The 
second  term  is  related  to  thje  electrostatic  interaction  between  the  emoedded 
atoms,  and  the  third  term  is  the  "band  structure"  energy  related  to  the  change  in 
one-electron  eigenvalues  in  going  from  the  reference  system  to  the  real  one  The 
third  term  is  usually  small  in  magnitude  in  metallic  systems,  but  is  important  in 
determining  eg  the  trends  over  transition  metal  series  of  the  s-d  hybridisation 
terms.  The  decomposition  (l)  bears  resemblance  to  the  familiar 
Dseuoopotential-based  Derturbation  theory for  simple  rnetais,  but  is  mucn 


wider  in  applicability.  There  are  several  formulations  based  on  ( 1  )  wmcn  differ 


in  detail  but  conform  to  the  same  philosophy. 

There  is  an  ongoing  effort  to  calculate  the  individual  terms  in  the 
expansion  (1  )  from  first  principles  -  with  a  minimum  of  adjustable  parametrs 
and  uncontrolled  approximations.  In  the  formulation  of  Jaccosen  et  j/16,  the 
function  F  is  based  on  the  homogeneous  electron  gas  and  reads 

Fin)  =  AE!',om(n)  -o<n,  (3) 

wnere 

-  - 1  dr  A^(r) .  (4) 

Above,  AEhom  is  the  embedding  energy1 '  of  an  atom  into  a  homogeneous  electron 
gas  of  density  n,  and  a <f>  is  the  atom-induced  change  in  the  electrostatic 
potential  of  the  jellium.  inserting  (3)  into  Eq.  (I)  implies  a  self-consistency 

requirement  in  constructing  the  pseudoatom  density  clouds:  na(R)  should  be 

obtained  from  embedding  the  atom  into  a  jellium  characterised  by  the  "tail 

summation”  density  nj(R)  By  going  to  the  atomic-sphere-approxirnation, 

•Jacobsen  et  a!  16  discuss  in  detail  the  various  cancellations  and  tail  corrections 
in  Eq.  (1).  As  emphasized  by  these  authors,  the  variational  property  of 
density-functional  theory  guarantees  that  independent  approximations  made  in 
both  the  electron  density  and  effective  potential  bring  about  only  second-order 


errors  in  total  energy.  This  property  is  useful  in  seeking  for  convenient 
approximations  in  the  evaluation  of  Eq.  ( 1 ). 

For  perfect  close-packed  monoatomic  solid  the  second  term  in  Eq.  ( I ) 
is  nearly  zero,  and  the  cohesive  energy  per  atom  reduces  to  a  very  simple 
formula 


Econ  =  AEhom(n)  -o<n  -  Ebs 
=  Ec(n)  ♦  EDs 

For  simple  metals  even  the  oand-structure  term  Eb5  is  small.  Thus  the  position 
and  depth  in  the  minima  of  Ec(n)  (see  Fig.  1.)  determine  the  equilibrium  lattice 
constant  and  cohesive  energy.  The  bulk  modulus  is 


(5) 

(6) 


B  =  (12  ttS)'1  d2Ec/d52  , 


where  5  is  the  Wigner-5eltz  radius.  The  equilibrium  configuration  is  in  general  a 
result  of  the  competition  between  inter-atomic  electrostatic  attraction  (  -ctfn) 
and  kinetic  energy  repulsion,  which  dominates  AEbom(n), 

Jacobsen  et  j/'°  have  also  shown  that  Eq.  (1),  with  appropriate 
account  of  the  electrostatic  second  term,  which  is  nonzero  in  a  low-symmetry 
situation,  can  easily  be  used  to  analyze  trends  and  even  quantitatively  predict 


energies  associated  with  ,  say,  surfaces  and  adsorbates. 


For  example, 


effective-medium  theory  offers  a  conceptually  simple  way  of  rationalizing  the 
ODservations  of  surface  relaxations  on  open  fee  (110)  surfaces  of  metals.  The 
driving  force  behind  the  inwards  relaxation  of  the  first  lattice  plane  is  the 
attempt  of  these  atoms  to  place  themselves  in  the  optimum  electron  density  as 

defined  bv  the  minimum  of  the  Er(n)  curve  of  Fig.  I.  The  atoms  in  the  first  plane 

have  fewer  nearest  neighbors  than  m  bulk  and  subsequently  sample  too  low  a 
density,  which  is  compensated  for  by  moving  closer  to  the  second  plane.  The 
electrostatic  term  (beyond  the  atomic-sphere-approximation)  opposes  this  trend 
and  diminishes  the  relaxation,  in  fact  this  term,  in  combination  with  the 
increased  density  from  the  relaxed  first  layer,  is  responsible  for  the  outwards 
relaxation  of  the  second  layer.  For  close-packed  surfaces  these  two  terms 
nearly  cancel,  but  the  more  open  the  surface  the  stronger  the  driving  force  from 

the  cohesive  term  Ec(n)  Similar  considerations  can  be  used  to  explain  the 

stability  of  the  missing-row  reconstruction  of  some  fee  (110)  surfaces. 
Likewise,  the  tendency  of  adsorbed  dimers  on  metals  of  (i)  to  have  a  shorter  bond 
length  than  in  vacuum  and  (ii)  to  rise  slightly  above  the  surface  as  compared  to 
single  adsorbates  can  be  easily  be  qualitatively  understood  as  arising  from  the 
same  basic  competition  between  the  embedding  cohesion  and  electrostatic 
repulsion. 

Redfield  and  Zanqwill'8  have  recently  used  the  information  contained 
;n  Fig.:  to  discuss  the  stability  of  icosahedrai  phases  in  Al-transition  metal 
alloys  Fig.  I  snows  that  Al  and  Mn,  for  example,  have  very  different  optimum 


density  requirements.  Redfield  and  Zangwill  argue  that  the  Mackay  icosanedron 
encased  in  a  sheil  of  A!  atoms  is  the  optimum  compromise  for  the  structure  of 
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the  alloy.  It  is  also  easy  to  see  from  Fig.  1 

that  icosanedral  binary  alloys  of  A!  are  only  formed  with  transition  metais  to 
the  right  of  Ti  in  the  Periodic  Table. 

2.3.  Embedded-atom  and  related  methods 

The  same  basic  physical  considerations  as  discussed  above  are 
involved  in  the  recent  ''embedded-atom"19  and  "glue"^  models  for  interatomic 
interactions.  The  total  energy  expression  in  these  is  very  similar  to  that  of  Eq. 
(1): 


Etot(N)=  IF 

However,  m  these  approaches  the  functions  Fj  and  <p  are  empirically 

determined,  obtained  by  choosing  reasonable  forms  for  the  two  functions  and 

fitting  to  experiment.  In  particular,  Daw  cYj/1  9  assumed  the  pair  potential  term 

<p  to  be  purely  repulsive 

?AB(R)  =  ^A^Lie^El  , 

R 


j(n j(R j ))  *£  2  <p(R|  -Rj).  (8) 


(9) 


where  ZA  and  are  vQ  -decendent)  effective  charges  i*  0)  of  two  err.Oeddod 


atoms  A  and  B  separated  Dy  distance  R.  The  "many-Dody"  potential  function  F(n), 
dependent  again  on  the  electron  density  atom  the  atom  positions,  starts  from 
zero  v/ith  a  negative  slope  and  positive  curvature.  It  can  be  obtained  by  fitting'-' 
the  expression  (  3)  to  give  cne  correct  sublimation  energy,  lattice  constant, 
elastic  constants,  vacancy  formation  energy  etc.,  whereafter  the  model  C3n  be 
used  m  computer  simulations  of  other  more  complicated  phenomena2-  The 
information  used  in  the  empirical  fit  actually  only  determines  F(n)  and  its  first 
two  derivatives  for  electron  densities  near  the  average  host  electron  density 

neq  Foiles42  has  arguea  that  the  zero  pressure  equations  of  state  can  be  used  to 

determine  F(n)  for  a  larger  class  of  electron  densities.  This  is  based  on  the 
ansatz  of  a  "universal"  equation  of  state  suggested  by  Rose  er  j/24  it  should  be 
noted,  however,  that  the  functions  F(n)  in  Eqs.  (I)  and  (3)  are  not  physically 
comparable,  the  latter  is  merely  one  possible  parametrisation  of  the  many-atom 
effects,  while  the  first  one  actually  obtains  the  function  from  embedding  into 
the  reference  system. 

Daw,  Foiles  and  coworkers2 ' ,22  have  applied  their  scheme  to  a  large 
variety  of  important  materials  science  problems.  The  computational  cost  does 
not  significantly  exceed  that  of  pair  potential  simulations,  once  the  functions  F 
and  <f>  are  fitted.  Ercolessi  et  j/.  20  have  used  a  very  similar  "glue"  model  to 
mvestigate  the  surface  reconstruction  of  Au(  100) 

Finms  and  Sinclair-4  have  suggested  yet  another  energy  expression 


whicn  resemblance  to  Eq.  ( I ).  in  particular,  they  choose 

F(n)  =  -  A  V7T  (10) 

Moreover,  n  in  their  work  is  not  interpreteo  as  a  density,  hut  the  rather  the 

superposition  is  thought  of  in  a  more  general  sense  that  na  in  Eg.  (2)  is  an 

adjustable  function.  The  square  root  in  Eq.  (10)  is  chosen  to  mimic  the  result  of 

tight-Dinding  theory,  in  which  n3  would  he  interpreted  as  a  sum  of  squares  of 

overlap  integrals.  This  recipe  and  fit  lead  to  an  N-Dody  potential  applicable  to 
bcc  transition  metals. 

3.  Computer  simulations  of  structural  properties 

The  above  ideas  have  recently  been  applied  in  a  number  of  computer 
simulation  studies  for  materials  problems  of  great  practical  interest.  These 
include  (i)  sputtering  and  associated  defect  production25  by  hyperthermal  ions 
near  single-crystal  surfaces,  (ii)  quenched-in  structures26  and  collision 
cascades  in  amorphous  solids,  and  effects  associated  with  melting  and 
solidification  in  constrained  geometries  and  near  surfaces  and  interfaces-'  : 
choose  to  discuss  one  example  of  the  latter  the  in  some  detail  below 
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3.1.  Rare  gas  -  meta!  interfaces 

Noble  gases  are  practically  >nsoiuble  in  metals  and  therefore  snow  a 
strong  tendency  to  precipitate  into  wnat  are  conventionally  called  bubbles  even 
in  the  cases  where  tne  precipitate  is  in  a  solid  phase  The  pressure  inside  these 
cubbies  can  oe  enormous,  up  to  several  GPa.  The  precipitation  tendency  naturally 
arises  from  the  strongly  repulsive  nature  of  the  metal-  rare  gas  interaction  m 
order  to  correctly  characterize  the  gas  densities  and  pressures  inside  rnese 
bubbles  via  a  numoer  of  experimental  probes*-3  it  has  turned  out  to  be  important 
to  have  an  accurate  structural  description  of  the  interface  between  the  metal 
and  the  gas. 

we  have  undertaken  a  series  of  molecular  dynamics  simulations*^  T'or 
confined  rare-gas  bubbles,  in  general,  large  precipitates  have  three-dimensional 
faceted  shapes.  The  problem  of  describing  the  atomic  structure  can  therefore  oe 
approached  by  considering  planar  metal-noble  gas  interfaces  representing  the 
facets  A  sandwich  geometry  was  set  up  with  six  layers  of  metal  atoms  (up  to 
432  atoms;  on  each  side  of  a  collection  of  up  to  600  gas  atoms  Periodic 
boundary  conditions  were  applied  in  all  three  directions  The  interatomic 
potentials  were  specified  according  to  effective-medium  theory  in  addition,  a 
long-range  van  der  waals  attraction  adjusted  to  give  correct  phvsisorption 
properties  of  single  gas  atoms  was  added  to  the  force  law  The  met3l-metai 


potential  reproduces  yve!  1  tre  bulk  lattice  constant,  ouik  modulus  ana 


sublimation  energy.  The  simulations  were  run  up  to  1500  hundred  timesteps  of  5 


*i0-'5  s  after  the  system  was  first  equilibrated  to  the  required  temperature 

Fig.  2  shows  two  examples  of  the  perpendicular  gas  density  profiles 
for  He  near  3n  Al  ( 100)  interface  No  solidification  is  observed  in  this  case  with 

increasing  mean  density  nHe  defined  as  the  average  value  for  large  a  strong 

peak  in  the  He  density  profile  first  develops  near  the  Al  surface  and  secondary 
peaks  in  the  fluid  follow  The  surface  peak  originates  from  the  statistical 
packing  of  atoms  near  the  repulsive  metal,  much  as  with  hard  spheres  near  an 
impenetrable  wall  The  other  pe3ks  at  higher  fluid  densities  signal  incipient 
ordering  in  the  compressed  fluid.  The  distance  between  the  peaks  corresponds  to 
that  between  close-packed  planes.  No  appreciable  enhancement  of  the  'ater3i 
ordering  is  observed  in  the  He  layer  closest  to  the  interface  compared  to  tne 
deeper  layers 

Another  set  of  simulations  was  performed  for  the  Cu/Kr  interface  ;h 
this  case  a  solidification  transition  is  observed.  At  high  enough  densities,  tr,ere 
are  solid  structures  with  lattice  planes  parallel  to  but 

incommensurate  with  the  metal  planes  These  can  be  obtained  frcm  var,ous 
starting  arrangements  with  different  packing  densities  At  lower  gas  dens’ties, 
a  fluid-like  profile  develops  at  the  interface  The  observed  general  features  </ 
the  pnase  diagram  are  consistent  with  bulk  Kr  thermodynamical  data 

Perhaps  a  more  interesting  observation  concerns  an  "anneal mg'  ~jn. 
performed  at  Kr  density  of  25  *  10*2  cm'-3  The  fluid  was  f'rst  equilibrated  at 
a  high  temperature,  whereafter  temperature  was  lowered 


oelow  trie  solidification  line.  Time  evolution  of  the  system  was  then  followed  at 
a  fixed  temperature  In  this  case  solid  atomic  planes  nucleate  from  the  interface 
(Fig  3)  The  planes  of  Kr  which  start  to  condense  parallel  to  the  Cu  surfaces  are 
those  of  close-packed  (lit)  character  They  are  not,  however,  perfect  (111) 
planes,  presumaoly  due  to  the  constraint  of  the  periodic  Coundarv  conditions 

Tne  oDservaticn  that  the  close-packed  din  planes  of  the  gas  anon 
themselves  with  the  metal  planes  is  also  consistent  with  the  experimental 
results  of  electron  diffraction,  it  seems  that  the  maior  criterion  for  tne 
orientational  relationship  of  the  rare  gas  Duobles  is  that  the  close-packed  planes 
of  the  gas  are  parallel  to  those  of  the  host  matrix  (tne  latter  presumaoly  form 
the  facets  of  the  metal  voids) 

This  result  seems  at  first  quite  surprising,  since  the  gas-metal 
interactions  at  the  relevant  distances  are  highly  repulsive  One  mignt  nave 
expected  the  gas  to  precipitate  with  a  loosely  packed  plane  next  to  a 
close-packed  metal  surface  in  addition,  since  the  gas-gas  interactions  are  aiso 
repulsive  at  small  distances,  the  close-packed  gas  has  a  higher  surface  energy 
than  a  more  loosely-packed  one  Finms^0  has  shown  that  the  solution  to  the 
apoarent  paradox  lies  in  the  fact  that  a  lowering  of  the  bulk  internal  energy  can 
be  achieved  with  the  epitaxial  orientation,  which  more  than  outweighs  the  higher 
interface!  energy  under  high  gas  pressure  and  constant  volume  conditions  tne 
placing  of  low-mdex  fee  <  1 1 1 )  surfaces  next  to  tne  metal  facets  is  more  than 
compensated  by  tne  relaxation  of  the  bulk  gas  which  this  permits  While  this 
result  is  seemingly  'naeoenoent  of  tne  actual  orientation  of  the  metal  facets,  tne 


1% 


edge  effects  inherent  in  the  geometry  used  in  the  simulation  might  be  important 
m  establishing  the  rather  small  energy  differences  between  hep  and  fee  ordering. 

Finally,  there  is  an  interesting  connection  between  the  molecular 
dynamics  results  and  the  experimentally  observed  superheating  of  rare  gas 
precipitates  in  metals.  Roussow  and  Donnelly^1  reported  a  solid-fluid  transition 
in  Ar  oubbies  in  Al  480  K  above  the  Dulk  Ar  melting  temperature  for  the  density 
question  The  melting  transition  was  deduced  from  the  disappearance  of  electron 
diffraction  spots  associated  with  Ar.  (The  superheating  effect  is  apparently 
much  smaller  or  even  absent  for  Kr  in  Cu  and  Ni.)  Since  the  bubbles  in  general 
must  oe  small  to  withstand  the  high  pressures  needed  to  solidify  the  fluid  at 
high  temperatures,  a  large  fraction  of  the  atoms  will  be  at  the  bubble  surface. 
Roussow  and  Donnelly^1  discuss  the  surface  effects  in  terms  of  reduced  thermal 
vibrations  of  the  gas  atoms  caused  by  the  metal  surface  ,  which  they  present  as 
a  mecnamsm  for  suppression  of  the  melting  transition.  The  present  molecular 
dynamics  results  suggest  an  alternative  explanation  ,  which  may  not  described 
as  true  superheating  effect  (it  does  not  exclude  that  in  the  bubble  interior 
melting  takes  place  at  the  predicted  bulk  melting  temperature).  The  peaks  in 
tne  density  profile  even  in  the  fluid  phase  are  signature  of  ordering  processes  at 
the  interface,  which  in  a  small  bubble  might  yield  a  sol  id- like  diffraction 
pattern  above  the  bulk  melting  point  of  the  solid.  The  inverse  mechanism, 
surface-initiated  melting  or  oremeltmg  has  been  observed  in  a  number  of 
computer  simulations^  and  ,  more  recently,  also  in  ion  scattering 
experiments^  The  picture  that  is  emerging  from  these  studies  is  that  melting 
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is  initiated  by  a  continuous  lattice  instability  at  the  free  surface  and  that  the 
melting  proceeds  along  directions  of  high  packing  density.  This  seems  to  De  the 
case  for  a  wide  variety  of  interatomic  potentials,  including  the  angular 
three-body  ones  necessary  to  stabilise  the  tetrahedral ly  coordinated  structures 
of  Si  and  Ge-'4 


A.  Conclusions 


To  sustain  the  steady  progress  in  expanding  the  scope  and  accuracy  of 
electronic  structure  calculations  ,  improvements  are  needed  in  algorithms  used 
to  solve  the  coupled  problems  of  electronic  (quantum)  and  ionic  (classical) 
degrees  of  freedom.  The  recent  ideas  based  on  simulated  annealing  and  related 
methods  are  most  promising  in  this  respect,  and  may  well  revolutionize  and 
unify  the  traditionally  somewhat  separate  fields  of  electronic  structure 
calculations  and  classical  molecular  dynamics  simulations.  Also  the  questions  of 
comDleteness  and  manageability  of  basis  sets  in  large  calculations  are  subject 
to  new  scrutiny.  On  the  other  hand,  the  development  of  simple  but  transparent 
schemes  such  as  the  effective-medium  theory  and  related  approaches  is  most 
helpful  from  both  pragmatic  and  conceptual  points  of  view. 

Eventually  the  limiting  obstacle  in  detailed  understanding  of  chemical 
and  structural  problems  lies  in  the  level  of  description  of  electronic  exchange 
and  correlation  effects  While  the  local-density  approximation  of 
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density-functional  theory  has  proven  surprisingly  robust  and  reliable,  its 
validity  is  fundamentally  limited^  Efforts  to  improve  and  generalize  it  3re 

thus  very  worthwhile. 
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Figure  Captions 

Fig. i.  The  energy  Ec(n)  (see  Eq.  (6)  )  as  a  function  of  electron  density  n  for 
various  elements.  After  refs.  1 6  and  17. 

Fig.2.  The  simulated  perpendicular  density  profiles  for  pressurized  He  fluid  near 
Al(  1 10)  wall.  T  »  300  K  After  ref.  29. 


Fig.  3.  The  simulated  nucleation  as  a  function  of  time  (0  -  c)  of  solid  Krd  1 1) 
layers  near  a  Cud  00)  interface  after  quenching  from  the  melt  (a).  After  ref  29 
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